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ABSTRACT 

A  study  for  validating  a  time-accurate  explicit  finite 
element  code  for  modeling  fully-coupled  flexible 
multibody  systems  carrying  liquid-filled  tanks  is 
presented.  The  multibody  system  includes  rigid  bodies, 
flexible  bodies,  joints,  and  actuators.  Rigid  bodies 
rotational  equations  of  motion  are  written  in  a  body-fixed 
frame  with  the  total  rigid  body  rotation  matrix  updated 
each  time  step  using  incremental  rotations.  Flexible 
bodies  are  modeled  using  total-Lagrangian  spring,  truss, 
beam  and  hexahedral  solid  elements.  A  penalty  model  is 
used  to  impose  the  joint/contact  constraints.  An  asperity- 
based  friction  model  is  used  to  model  joint/contact 
friction.  The  fluid  governing  equations  of  motion  are  the 
incompressible  Arbitrary  Lagrangian-Eulerian  Navier- 
Stokes  equations  along  with  a  large-eddy  simulation 
(LES)  turbulence  model.  The  fluid's  free-surface  is 
modeled  using  an  acceptor-donor  volume-of-fluid  based 
algorithm.  Coupling  between  the  fluid  and  solid  is 
achieved  by  solving  Newton’s  equations  of  motions  at 
the  fluid-solid  interface  nodes. 

The  validation  study  is  conducted  using  a  multibody 
system  consisting  of  a  rigid  baffled  tank  mounted  on 
suspension  springs.  The  springs  are  connected  to  a  rigid 
frame  mounted  on  two  linear  hydraulic-actuators. 
Experiments  with  various  input  ramp  and  harmonic 
excitation  from  the  actuators  are  performed  and  the 
results  of  the  experiments  are  compared  to  the  results 
obtained  using  the  model.  The  system  response  is 
measured  using  linear-displacement  transducers  at  the 
springs  and  two  cameras  showing  side  and  front  views 
of  the  tank.  The  results  show  that  the  model  can  predict 
with  reasonably  good  accuracy  the  test  system’s 
dynamic  response. 

1.  INTRODUCTION 

Many  practical  applications  involve  a  flexible  multibody 
system  carrying  one  or  more  liquid  filled  tanks.  The 
multibody  system  can  be  a  ground  vehicle  (truck,  train, 
or  car),  ship,  airplane  (commercial  jet,  military  jet,  or 
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helicopter),  a  space  vehicle  (rocket  or  reusable  launch 
vehicle)  or  a  space  structure  (space  station  or  satellite). 
The  tank  can  be  a  payload  tank  or  a  liquid  fuel  tank.  In 
those  applications,  an  accurate  computational  model  for 
predicting  the  coupled  solid-fluid  system  response  can 
greatly  reduce  the  cost  and  time  required  to  reach  an 
optimum  system  design  that  satisfies  the  various 
operating  constraints.  The  model  must  accurately 
account  for  the  following  effects: 

•  Incompressible  liquid  flow  in  a  moving/deforming 
container. 

•  Modeling  of  the  liquid  free-surface. 

•  Modeling  of  turbulence  and  viscous  effects. 

•  Coupling  between  the  solid  and  the  fluid  at  the  fluid- 
structure  interface. 

•  Large  rotation  of  the  solid  bodies. 

•  Deformation  of  the  solid  flexible  bodies.  Flexible 
bodies  can  be  modeled  as  beam,  shells  or  general 
solids.  For  example  a  shell  model  can  be  used  for  a 
flexible  tank. 

•  Joints  kinematic  constraints  including  joint  friction 
and  clearances. 

•  Frictional  contact.  For  example  for  ground  vehicle 
applications  the  rolling  frictional  contact  of  tires  need 
to  be  accurately  modeled. 

•  Actuators  and  control  laws. 

•  Motion  control  components  including  transmission 
components,  clutches,  and  brakes  for  ground 
vehicle  applications.  All  those  components  involve 
friction. 

The  fluid  flow  is  governed  by  the  incompressible  Navier- 
Stokes  equations.  Finite  volume  [1-3],  finite  element  [4- 
6]  or  particle  [7]  discretization  techniques  have  been 
used  to  model  the  fluid  flow.  Many  techniques  for 
modeling  fluid  flow  with  a  free-surface  have  been 
developed  in  the  literature.  They  include: 

(a)  Volume-of-fluid  (VOF)  method  [8,  9],  Each  element 
has  a  VOF  value  between  0  (for  empty  elements) 
and  1  (for  elements  completely  filled  with  fluid).  The 
free  surface  is  reconstructed  for  each  element  using 
piecewise-linear  planar  segments  that  are  calculated 
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from  the  VOF  value  of  the  element  along  with  the 
VOF  values  of  neighboring  elements  (which  are 
used  to  determine  the  normal  to  the  planar  surface). 
When  VOF  algorithms  were  first  used  the  free 
surface  was  reconstructed  using  either  vertical  or 
horizontal  surfaces  [8,  3],  The  VOF  values  of  all  the 
elements  are  updated  each  time  step  by  calculating 
the  mass  flux  between  elements.  The  mass  flux  for 
free-surface  elements  is  calculated  by  taking  into 
account  the  smaller  surface  through  which  the  fluid 
can  move  due  to  the  presence  of  the  free  surface. 

(b)  Level-set  method.  This  method  uses  a  smooth 
scalar  function  defined  at  every  node  in  the  fluid 
domain  which  specifies  the  signed  smallest 
Eucledian  distance  between  the  node  and  the 
interface  [5],  The  evolution  of  the  scalar  function  is 
governed  by  a  convection  transport  equation  where 
the  interface  is  moved  with  the  fluid  velocity. 

(c)  Arbitrary  Lagrangian-Eulriari  (ALE)  method.  Using 
this  method,  the  fluid  mesh  deforms  and  moves 
along  with  the  fluid’s  free-surface  [4,  6,  10],  A 
disadvantage  of  this  method  is  that  it  does  not  allow 
large  surface  deformation  including  surface  break¬ 
up  and  merging  unless  the  fluid  domain  is  re¬ 
meshed  [11,  12],  If  re-meshing  is  used  frequently,  it 
can  degrade  the  solution  accuracy  due  to  re¬ 
interpolation  of  the  solution  field  onto  the  new  mesh. 

(d)  Lagrangian  particle  methods.  Lagrangian  particles 
which  represent  small  packets  of  fluid  are  used  to 
model  the  fluid  flow.  A  contact  model  between  the 
fluid  particles  is  used  to  model  the  fluid 
compressibility  and  viscous  effects.  This  method 
naturally  handles  free-surfaces.  The  main 
disadvantage  of  those  types  of  methods  is  the  large 
number  of  particles  and  computing  time  needed  to 
accurately  predict  the  fluid  motion.  A  special  type  of 
this  class  of  methods  which  has  been  successfully 
applied  to  free-surface  flows  and  fluid-structure 
interaction  problems  is  the  particle  finite  element 
method  (PFEM)  [7]  in  which  the  particles  are  used  to 
generate  a  polyhedral  finite  element  mesh  every 
time  step  using  an  extended  Delaunay  tesselation. 
The  solution  of  the  incompressible  Navier-Stokes 
equations  is  then  carried  using  that  mesh.  The 
PFEM  requires  less  fluid  particles,  however,  the 
tessellation  step  is  computationally  intensive. 
Particle  methods  can  easily  handle  surface  break-up 
and  merging,  floating  bodies,  and  fluid-solid 
impenetrability  boundary  condition. 

Techniques  to  handle  fluid  flow  in  a  moving/deforming 

container  include: 

1.  Fixed  Cartesian  fluid  mesh  with  cut-cell  boundary 
condition.  The  fluid  domain  is  a  Cartesian  mesh.  The 
container  moves  inside  this  mesh.  A  cut-cell 
technique  is  used  to  find  where  the  boundary  of  the 
container  intersects  with  the  fluid  cells.  The  cut-cell 
surfaces  are  then  used  to  impose  wall 
impenetrability  and  adhesion  boundary  conditions  [1, 
2],  The  fixed  Cartesian  mesh  has  the  advantage  that 


mesh  generation  is  straightforward.  In  addition,  it 
allows  modeling  floating  objects  and  tank  baffles 
with  no  additional  effort.  However,  the  main 
disadvantage  of  the  technique  is  that  the  fluid-solid 
impenetrability  and  no-slip  boundary  conditions  are 
satisfied  only  in  a  time  average  sense.  Also,  the 
method  has  stability  and  accuracy  problems  when 
the  cut-cell  elements  at  the  solid-fluid  interface 
become  small. 

2.  Moving  ALE  mesh.  The  fluid  is  modeled  using  a  fluid 
mesh  that  moves  and  deforms  with  the  tank. 

3.  Fixed-fluid  mesh  with  the  Navier-Stokes  equations 
written  in  a  reference  frame  fixed  to  the  tank  [4,  13], 
Since  the  tank  frame  is  a  non-inertial  frame 
(accelerating  frame),  writing  the  equations  of  motion 
with  respect  to  that  frame  results  in  a  complex  inertia 
operator  which  involves  centrifugal  and  coriolis 
acceleration  terms.  Also,  this  method  cannot  -  by 
itself  -  deal  with  a  deforming  container. 

4.  Particle  methods.  Normal  contact  constraint 
between  the  particles  and  the  tank  wall  is  used  to 
model  the  wall  impenetrability  constraint.  Tangential 
friction  between  the  particles  and  the  tank  wall  is 
used  to  model  wall  adhesion  and  viscous  effects. 

In  order  to  model  fluid  flow  with  a  free-surface  in  a 
moving  deforming  container,  the  above  methods  must 
be  combined.  Table  1  shows  the  references  where  the 
various  combinations  of  the  above  techniques  were 
used. 

Table  1  References  for  the  various  combinations  of  techniques  for 
modeling  a  free-surface  flow  in  moving  deforming  container. 


In  the  present  paper,  an  ALE  mesh  is  used  for  modeling 
the  moving/deforming  container  along  with  a  VOF  free- 
surface  model.  This  method  does  not  suffer  from  the  cut¬ 
cell  solid-fluid  interface  boundary-condition  problem  like 
the  fixed  grids  with  cut-cell  boundary  condition  methods. 
It  also  allows  modeling  surface  break-up  and  merging 
without  the  need  for  re-meshing. 

A  review  of  multibody  dynamics  modeling  techniques 
including  deformation  reference  frames,  treatment  of 
large  rotations,  discretization  techniques,  finite  elements, 
constraint  and  contact  modeling,  and  solution 
techniques  is  presented  in  [14],  In  the  present  paper,  a 
flexible  multibody  dynamics  code  with  the  following 
characteristics  is  used: 


An  explicit  time-integration  solver  accurate  for  long 
simulation  times  [15]. 

Total  Lagrangian,  total  displacement  equations  of 
motion  formulation  with  the  degrees  of  freedom 
referred  to  a  global  inertial  reference  frame  [15-18]. 

A  library  of  truss,  beam,  and  solid  nonlinear  finite 
elements  with  Cartesian  coordinate  degrees  of 
freedom  allowing  arbitrarily  large  element  rotations. 
Those  include: 

o  Torsional-spring  type  3-node  beam 

elements  [15,  16]. 

o  Natural-modes  eight-node  brick  elements 

[17,  18].  Those  elements  can  also  be  used 
to  model  shells  and  beams.  One  element 
through  the  thickness  is  sufficient  to 
accurately  model  the  membrane,  shear,  and 
bending  characteristics.  They  do  not  exhibit 
locking  or  spurious  modes  (widely  used 
techniques  to  alleviate  locking  such  as 
hourglass  control  lead  to  elements  that  do 
not  maintain  solution  accuracy  over  very 
long  solution  times).  They  are 
computationally  efficient.  Assumed  strain 
elements  of  comparable  accuracy  are  more 
computationally  expensive.  Any  material  law 
can  be  used  with  those  elements  including: 
linear  elastic,  hyper-elastic,  and  non-linear 
laws. 

Penalty  formulation  for  modeling  joints  (spherical, 
revolute,  cylindrical  and  prismatic)  [19]. 

Rigid  body  rotational  equations  of  motion  are  written 
in  a  body  (material)  frame,  with  the  resulting 
incremental  rotations  added  to  the  total  body  rotation 
matrix  [19]. 

Normal  contact  modeled  using  a  penalty  formulation 

[20,  21]. 

Frictional  contact  modeled  using  an  accurate  and 
efficient  asperity-based  friction  model  [22], 

Special  elements  for  modeling  wheels/pulleys  [20], 
sprockets  [23],  and  clutches  [24], 

General  tire  model  [21].  The  model  includes  the 
details  of  the  tire  construction.  The  tire  rubber  is 
modeled  using  brick  elements.  The  bead,  tread  and 
ply  are  modeled  using  beam  elements  along  the  tire 
circumference  and  meridian  directions  with 
appropriate  stiffness  and  damping  properties. 
Normal  contact  between  the  tire  and  the  wheel  and 
between  the  tire  and  the  pavement  is  modeled  using 
the  penalty  technique.  Friction  is  modeled  using  an 
asperity-based  approximate  Coulomb  friction  model. 
The  tire  inflation  pressure  is  modeled  by  applying  a 
force  normal  to  the  inner  surface  elements  of  the  tire 
with  the  out-of-equilibrium  force  and  moment  applied 
to  the  wheel  to  guarantee  self-equilibrium  of  the  tire 
and  wheel  under  the  pressure  load  [21]. 


•  General  contact  search  algorithm  that  finds  the 
contact  penetration  between  finite  elements  and 
other  elements  as  well  as  general  triangle  and 
quadrilateral  surfaces. 

Two-way  coupling  between  the  multibody  system 
(vehicle)  motion  and  the  fluid  is  achieved  by  satisfying 
the  following  conditions  at  the  solid-fluid  interface: 

•  The  fluid  velocity  normal  to  the  solid’s  surface  must 
be  equal  to  the  normal  solid  velocity. 

•  The  fluid  velocity  tangent  to  the  solid  surface  can 
range  from  being  equal  to  the  tangential  velocity  of 
the  solid  surface  (no  slip  condition)  to  being  free. 

•  No  additional  energy  or  momentum  to  the  system 
should  be  introduced  at  the  interface. 

The  above  boundary  conditions  are  satisfied  by  using 
Newton’s  equations  of  motion  to  find  a  common  normal 
acceleration  for  the  fluid  and  the  solid  at  the  interface. 
The  tangential  fluid  and  solid  accelerations  can  range 
from  being  the  same  (no-slip  condition)  to  being 
completely  decoupled. 

In  the  present  paper,  a  single  computational  code  which 
uses  a  time-accurate  explicit  solution  procedure  is  used 
to  solve  both  the  solid  and  fluid  equations  of  motion. 
Many  commercial  software  and  studies  on  modeling 
liquid  sloshing  coupled  with  solid  body  motion  use  two 
codes  which  pass  the  interface  forces  and  motion  back 
and  forth  and  iterate  on  the  two  codes  until  equilibrium  is 
achieved  [e.g.  3,  13].  This  approach  adds  additional 
computational  burden  and,  in  general,  does  not  achieve 
the  same  accuracy  as  the  single  integrated  code 
solution  due  to  the  difficulty  in  achieving  an  equilibrium 
solution  between  two  disjoint  codes. 

The  finite  element  code  is  validated  by  comparing  its 
response  prediction  with  the  results  from  a  multibody 
system  consisting  of  a  rigid  baffled  tank  mounted  on 
suspension  springs.  The  springs  are  connected  to  a  rigid 
frame  which  is  mounted  on  two  linear  hydraulic- 
actuators.  Using  the  actuators,  the  frame  can  be  moved 
in  pitch,  roll  and  combination  pitch/roll.  A  test  matrix  with 
different  combinations  of  frame  motions,  tank  fill  levels 
(empty  or  half-full),  and  excitation  types  (ramp  or 
harmonic)  is  carried  out.  The  system  response  is 
measured  using  linear-displacement  transducers  at  the 
suspension  springs  and  two  cameras  showing  side  and 
front  views  of  the  tank  along  with  the  free-surface 
motion.  Comparisons  of  the  experimental  results  with 
the  results  generated  using  the  finite  element  code  are 
presented. 

The  rest  of  this  paper  is  organized  as  follows.  In 
Sections  2  and  3  the  equations  of  motion  for  the  solid 
and  fluid  are  presented.  In  Section  4  the  fluid-structure 
coupling  model  is  presented.  In  Section  5  the  VOF  free- 
surface  model  is  presented.  In  Section  6  the  overall 
explicit  solution  procedure  is  outlined.  In  Section  7  the 
validation  study  along  with  the  experimental  setup  are 
presented.  Finally,  in  Section  8  some  concluding 
remarks  are  offered. 


2.  SOLID  EQUATIONS  OF  MOTION 

In  the  subsequent  equations  the  following  conventions 
will  be  used: 

•  The  indicial  notation  is  used. 

•  The  Einstein  summation  convention  is  used  for 
repeated  subscript  indices  unless  otherwise  noted. 

•  Upper  case  subscript  indices  denote  node  numbers. 

•  Lower  case  subscript  indices  denote  vector 
component  number. 

•  The  superscript  denotes  time. 

•  A  superposed  dot  denotes  a  time  derivative. 

The  translational  equations  of  motion  are  written  with 
respect  to  the  global  inertial  reference  frame  and  are 
obtained  by  assembling  the  element  equations.  The 
finite  elements  used  here  use  only  translational  DOFs 
with  no  rotational  DOFs.  This  is  advantageous  in  terms 
of  computational  efficiency,  accuracy,  and  robustness 
[14].  Those  equation  also  include  the  rigid-body  (such  as 
the  wheel)  translational  DOFs.  The  equations  can  be 
written  as: 

Mkx'b=F.‘b .+Fa'B  (1) 

where  f  is  the  running  time,  K  is  the  global  node  number 
(no  summation  over  K;  K=1->/V  where  N  is  the  total 
number  of  nodes),  /  is  the  coordinate  number  (/=  1 ,2,3),  a 
superposed  dot  indicates  a  time  derivative,  MK  is  the 
lumped  mass  of  node  K,  x  is  the  vector  of  nodal 
Cartesian  coordinates  with  respect  to  the  global  inertial 
reference  frame,  and  x  is  the  vector  of  nodal 
accelerations  with  respect  to  the  global  inertial  reference 
frame,  Fs  is  the  vector  of  internal  structural  forces,  and 
Fa  is  the  vector  of  externally  applied  forces,  which 
include  surface  forces  and  body  forces. 

For  each  rigid  body,  a  body-fixed  material  frame  is 
defined.  The  rigid  body  is  represented  by  one  node 
located  at  the  body’s  center  of  mass,  which  is  also  the 
origin  of  this  frame.  The  mass  of  the  body  is 
concentrated  at  the  node  and  the  inertia  of  the  body 
given  by  the  inertia  tensor  Iy  is  defined  with  respect  to 
the  body  frame.  The  orientation  of  the  body-frame  is 
given  by  R%  which  is  the  rotation  matrix  relative  to  the 

global  inertial  frame  at  time  t0.  The  rotational  equations 
of  motions  are  written  for  each  rigid  body  with  respect  to 
its  body-fixed  material  frames  as: 

iJkj  =  T‘Ki  +  Ki  -  fe,  x  (IkAj ))*  (2) 

where  IK  is  the  inertia  tensor  of  rigid  body  K,  $  and  $ 

are  the  angular  acceleration  and  velocity  vectors’ 
components  for  rigid  body  K  relative  to  it’s  material 
frame  in  direction  j,  TsKi  is  the  component  of  the  vector  of 
internal  torque  at  node  K  in  direction  i,  and  TaKi  is  the 
component  of  the  vector  of  applied  torque  in  direction  i. 
The  summation  convention  is  used  only  for  the  lower 
case  indices  i  and  j. 


The  trapezoidal  rule  is  used  as  the  time  integration 
formula  for  solving  equations  (1)  for  the  global  nodal 
positions  x: 


Xxf*  +  0.5  At  (xlKj  +  x*KjAt ) 

(3a) 

xKjAt  +  0-5  At  (xfKj  +  xfKjAt) 

(3b) 

where  At  is  the  time  step.  The  trapezoidal  rule  is  also 
used  as  the  time  integration  formula  for  the  nodal 
rotation  increments: 

e‘Kj  =  e^l+o.5M{elKj  +  et-;t)  (4a) 

A6‘Kj  =  0.5  A?  {0‘Kj  +  e‘~^)  (4b) 

where  A6Kj  are  the  incremental  rotation  angles  around 
the  three  body  axes  for  body  K.  The  rotation  matrix  of 
body  K  (Rk)  is  then  evaluated  using: 

R^R^Ri  A0fKj)  (5) 

where  R( A0fKj)  is  the  rotation  matrix  corresponding  to 

the  incremental  rotation  angles  from  Equation  (4b). 

The  explicit  solution  procedure  used  for  solving 
equations  (1-5)  along  with  constraint  equations  is 
presented  in  Section  6.  The  constraint  equations  are 
generally  algebraic  equations,  which  describe  the 
position  or  velocity  of  some  of  the  nodes.  They  include: 

•  Prescribed  motion  constraints: 

f({x},t)  =  0  (6) 

•  Joint  constraints: 

/«*})  =  0  (7) 

•  Contact/impact  constraints: 

/(M)>0  (8) 

The  penalty  technique  is  used  for  imposing  the 
constraints  in  which  a  normal  reaction  force  is  generated 
when  a  node  penetrates  into  a  contact  body.  The 
magnitude  of  the  force  is  proportional  to  the  penetration 
distance  [20,  21].  An  asperity-spring  friction  model  is 
used  to  model  joint  and  contact  friction  [22]  in  which 
friction  is  modeled  using  a  piece-wise  linear  velocity- 
dependent  approximate  Coulomb  friction  element  in 
parallel  with  a  variable  anchor  point  spring.  The  model 
approximates  asperity  friction  where  friction  forces 
between  two  rough  surfaces  in  contact  arises  due  to  the 
interaction  of  the  surface  asperities. 

3.  SEMI-DISCRETE  FLUID  EQUATIONS  OF 
MOTION 

The  dynamic  response  of  the  fluid  is  described  by  the 
ALE  version  of  the  incompressible  Navier-Stokes 
equations,  namely,  the  equations  of  conservation  of 
momentum  and  mass  for  a  moving  deforming  control 
volume  along  with  a  large-eddy  simulation  (LES) 
transport  equation.  Those  are: 
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(15a) 

(15b) 


where  V  is  the  element  volume,  t  is  the  running  time,  p  is 
the  density  of  the  fluid,  u  is  the  fluid  velocity  vector 
relative  to  the  global  reference  frame,  *5  is  the  fluid 
velocity  vector  relative  to  the  moving  fluid  mesh,  v  is  the 
velocity  of  the  fluid  mesh,  P  is  the  relative  pressure,  x  is 
the  position  vector,  r  is  the  deviatoric  stress  tensor,  D  is 
the  rate  of  deformation  tensor,  /  is  the  body  force 

vector,  r  is  the  artificial  compressibility  parameter,  p0  is 
the  nominal  fluid  density,  //  is  the  fluid  viscosity,  /*  is  an 
additional  turbulence  viscosity  calculated  using  Equation 
15b),  K  is  the  eddy  kinetic  energy,  and  s  is  the  sub-grid 
scale  eddy  kinetic  energy  dissipation  term. 
Incompressible  flow  is  modeled  using  the  artificial 
compressibility  technique  [25],  A  finite  element 
formulation  is  used  to  derive  the  element’s  semi-discrete 
equations  of  motion  from  the  governing  equations  (9- 
13).  8-node  hexahedral  elements  are  used  with  tri-linear 
equal-order  velocity  and  pressure  interpolation.  A 
pressure  averaging  algorithm  [26]  is  used  to  eliminate 
pressure  checker-boarding  (due  to  the  use  of  an  equal 
order  interpolation  for  pressure  and  velocity).  The 
element  equations  are  assembled  into  the  global  semi¬ 
discrete  equations  of  motion: 


II 

-  ^ 

•s 

(16) 

VtN  P‘Ni  =  Q/n 

(17) 

MfKK‘N=Sf‘N 

(18) 

pkj  =  pkJa>  +  0.5  Ar  (P‘Kj  +  P‘k~m) 
K'kj  =  K‘~^  +  0.5  Ar  (k‘Kj  +  K‘~^ ) 


4.  FLUID-STRUCTURE  COUPLING  MODEL 

Newton’s  equations  of  motion  are  used  to  find  a 
common  normal  acceleration  for  the  fluid  and  the  solid  at 
the  interface.  This  is  done  for  each  node  at  the  fluid- 
structure  boundary  as  follows: 

(ms  +  m/)un  =  y  Fluid  Forces  +  y  StructureForces  (22a ) 

(ms + m/)vn  =  y  Fluid  Forces + y  StructureForces  (22b) 

where  ms  is  the  solid  mass  of  the  node,  mt  is  the  fluid 
mass  of  the  node,  ii,,  and  v„  are  respectively  the  fluid 
and  solid  accelerations  of  the  node  normal  to  the  fluid- 
structure  interface.  The  tangential  fluid  and  solid 
accelerations  ( {/, ,  i>, )  are  calculated  using  the  following 
equations: 

((1  - s)ms  +  mr)ut  =  (1- s)^StructureForces +^Fluid Forces  (23a) 

(ms  +  (1  -  s)mi)i>,  =  y  StructureForces  +  (1  -  s)^  Fluid  Forces  (23b) 

where  s  is  the  slip  factor.  A  no-slip  condition 
corresponds  to  a  slip  factor  of  zero.  The  slip  factor 
determines  how  much  of  the  fluid  and  structure  forces 
are  mutually  exchanged.  Equations  22  and  23  are 
written  for  all  fluid-solid  interface  nodes.  The  fluid  mesh 
must  move  and  deform  with  the  tank.  This  is  done  by 
modeling  the  fluid  mesh  using  very  light  and  compliant 
(3  orders  of  magnitude  less  than  the  tank)  solid  brick 
elements  (called  “mock”  mesh).  The  ALE  formulation  is 
used  to  account  for  the  fluid  mesh  deformation/motion. 

5.  VOF  FREE-SURFACE  MODEL 

For  each  fluid  element  a  VOF  value  between  0  and  1  is 
defined,  where  0  corresponds  to  empty  elements  and  1 
corresponds  to  elements  completely  filled  with  fluid.  The 
elements’  VOF  values  are  updated  each  time-step  by 
moving  fluid  from  a  completely  or  partially  filled  “donor” 
element  to  an  empty  or  partially  filled  neighboring 
“acceptor”  element  using  the  following  model: 


where  MtN  is  the  lumped  fluid  mass  of  node  N,  um  is 
component  /  of  the  fluid  acceleration  at  node  N,  Ff‘m  is 

component  /  of  the  fluid  forces  at  node  N,  VfN  is  the 
lumped  fluid  volume  at  node  N,  p‘N  is  the  fluid  pressure 

rate  at  node  N,  Qf'N  is  the  fluid  pressure  flux  at  node  N, 

K‘n  is  the  eddy  kinetic  energy  rate  at  node  N,  and  S/N  is 

the  eddy  kinetic  energy  flux  at  node  N.  Those  equations 
are  integrated  using  the  trapezoidal  rule  along  with  an 
explicit  solution  procedure  to  yield  the  nodal  fluid  velocity 
and  pressure: 

ukj  =  11  /,/V  +  9-5  At  (u‘Kj  +  u‘KjAt)  (1 9) 


Veo  =  Ve  VOF, 


Vna  =  Vn  (1~  VOFn) 

[A tSAn.u  Veo>  AV  and 


AV  =  < 


Veo 

Vna 


Veo  <  AV 
Vna  <  AV 


Vna  >  Veo 


(24) 

(25) 

(26) 


where  Ve  is  the  volume  of  the  element;  V„  is  the  volume 
of  the  neighboring  element;  Veo  is  the  volume  of  the 
element  occupied  by  the  fluid;  Vna  is  the  volume  of  the 
neighboring  element  available  to  receive  fluid;  AV  is  the 
volume  flow  through  the  boundary  between  the  two 
elements  in  a  time  step;  At  is  the  solution  time  step;  S  is 
the  surface  area  between  the  two  elements;  A  is  a  value 


between  0  and  1  indicating  the  free-surface  aperture 
through  which  the  fluid  can  move  from  the  element  to 
the  neighboring  element;  «  is  a  unit  vector  normal  to  S; 
and  u  is  the  fluid  velocity  vector  at  the  surface  S.  If  AV  is 
less  than  0  then  the  element  is  an  acceptor  element  and 
the  VOF  values  are  not  updated  because  they  will  be 
updated  later  when  the  neighbor  element  is  set  to  be  the 
donor  element.  If  AV  is  greater  than  0  then  the  VOF 
values  are  updated  using  the  following  equations: 

VOFe  =  VOFe  -AV/Ve  (27) 

VOFn  =  VOFn  +  A  V/Vn  (28) 

The  free-surface  apertures  A  at  the  element  interfaces 
are  used  to  limit  the  fluid  flow  based  on  the  location  of 
the  free  surface  inside  the  element.  A  is  calculated  as 
follows.  If  the  VOF  value  of  the  element  is  1  then  there  is 
no  free-surface  at  the  element,  therefore  A= 1.  For 
elements  with  a  VOF  value  less  than  1,  the  following 
steps  are  used  to  calculate  A: 

•  Calculate  the  normal  to  the  surface  by  looking  at  a 
stencil  of  neighboring  elements  around  the  element. 
This  is  done  using  the  following  equation: 

net  =  VOF„,  Snk  nnu  i=  1,2,3  (29) 

where  «e,  is  the  /th  component  of  the  normal  to  the 
free-surface  at  the  element,  VOFni  is  the  VOF  value 
for  neighboring  element  number  k,  Snk  is  the  area  of 
the  intersection  surface  between  the  element  and 
neighboring  element  k,  and  n„ki  is  the  component  i  of 
the  normal  to  the  surface  between  the  element  and 
neighboring  element  k.  Tie  is  then  normalized  into  a 
unit  vector.  Figure  1  shows  a  2D  4-node  quadrilateral 
and  the  free-surface  along  with  the  normal  Tu  . 

•  Calculate  the  apertures  A  for  each  neighboring 
element  by  constructing  a  planar  surface  with 
normal  Tie  and  with  total  volume  equal  to  VOFe  Ve 
(see  Figure  1). 


Figure  1  Stencil  of  neighboring  elements  used  to  determine  the  free- 
surface  normal  ne  and  the  liquid  free-surface. 

6.  EXPLICIT  SOLUTION  PROCEDURE 

The  solution  fields  for  modeling  the  solid,  fluid  and  liquid 
free-surface  are  defined  at  the  model  nodes.  These  are: 


•  Solid  translational  positions. 

•  Solid  translational  velocities. 

•  Solid  translational  accelerations. 

•  Solid  rotation  matrices. 

•  Solid  rotational  (angular)  velocities. 

•  Solid  rotational  (angular)  accelerations. 

•  Fluid  velocities. 

•  Fluid  accelerations. 

•  Fluid  pressure. 

•  Fluid  pressure  rate. 

•  Volume-of-fluid. 

•  Eddy  kinetic  energy 

•  Eddy  kinetic  energy  rate. 

The  explicit  time  integration  solution  procedure  for 
modeling  the  coupled  response  of  the  solid  (multibody 
system),  fluid,  and  liquid  free-surface  (using  the  VOF 
formulation)  predicts  the  time  evolution  of  the  above 
response  quantities.  The  procedure  is  implemented  in 
the  DIS  [27]  (Dynamic  Interactions  Simulator) 
commercial  software  code  and  is  outlined  below: 

1)  Prepare  the  run: 

a.  Set  the  initial  conditions  for  the  solution  fields 
identified  above. 

b.  Create  a  list  of  all  the  finite  elements  (including 
both  solid  and  fluid  elements). 

c.  Create  a  list  of  elements  that  will  run  on  each 
processor.  This  is  done  using  an  algorithm  which 
tries  to  make  the  computational  cost  on  each 
processor  equal. 

d.  Create  a  list  of  all  the  constraints  (including  both 
solid  and  fluid  constraints). 

e.  Calculate  the  solid  masses  for  each  finite  element 
node  by  looping  through  the  list  of  finite  elements. 
Note  that  the  solid  masses  are  fixed  in  time. 

f.  For  each  node  create  a  list  of  corner  and  edge 
nodes  that  are  connected  to  it  using  fluid  volume 
elements. 

g.  VOF  preparations: 

i.  Find  a  list  of  the  volume  fluid  elements. 

ii.  Create  a  list  of  fluid  volume  elements  that  will 
run  on  each  processor.  This  is  done  using  an 
algorithm  which  tries  to  make  the  computational 
cost  on  each  processor  equal. 

iii.  For  each  element  find  all  neighboring  elements. 

iv.  For  each  element  find  the  element  VOF  using 
the  nodal  VOF. 

v.  Re-interpolate  the  elements’  VOF  to  nodal  VOF. 

h.  Loop  over  all  the  elements  and  find  the  minimum 
time  step  for  the  explicit  solution  procedure. 

i.  Loop  over  all  the  elements  and  create  a  list  of  wall 
nodes.  For  each  wall  node  find  the  list  of  fluid 
boundary  elements. 

2)  Loop  over  the  solution  time  and  increment  the  time 
by  At  each  step  while  doing  the  following: 

a.  Set  the  nodal  values  at  the  last  time  step  to  be 
equal  to  the  current  nodal  values  for  all  solution 
fields. 

b.  Do  2  iterations  (a  predictor  iteration  and  a 
corrector  iteration)  of  the  following: 


i.  Initialize  the  nodal  fluxes  to  zero.  Those  include: 
solid  forces,  solid  moments,  fluid  forces, 
boundary  fluid  forces,  and  pressure  fluxes.  In 
addition,  the  lumped  nodal  fluid  volume  and  fluid 
mass  vectors  are  also  initialized  to  zero. 

ii.  Calculate  the  nodal  solid  and  fluid  fluxes  and  the 
lumped  fluid  volume/mass  vectors  by  looping 
through  all  the  elements  while  calculating  and 
assembling  the  element  nodal  fluxes  and 
vectors.  This  is  the  most  computational  intensive 
step.  This  step  is  done  in  parallel  by  running 
each  list  of  elements  identified  in  step  1  .c  on  one 
processor. 

iii.  Find  the  nodal  values  at  the  current  time  step 
using  the  semi-discrete  equations  of  motion  and 
the  trapezoidal  time  integration  rule  (Equations 
1-5  and  19-21). 

iv.  Execute  the  solid  and  fluid  constraints.  The 
constraints  prescribe  the  nodal  values. 

v.  Apply  fluid-structure  interface  boundary 
conditions  for  all  wall  nodes  found  in  Step  l.i 
(see  Equations  22,  23).  This  is  done  by  doing 
the  following  for  each  wall  node: 

1.  Find  the  normal  to  the  surface  at  the  wall 
node. 

2.  Normalize  the  surface  normal. 

3.  Find  the  solid,  fluid  bulk  and  fluid  boundary 
forces  in  the  directions  normal  and  tangent  to 
the  surface. 

4.  Find  the  normal  and  tangential  solid  and  fluid 
accelerations  using  the  trapezoidal 
integration  rule  and  the  wall  slip  percentage. 

vi. Set  the  pressure  boundary  conditions  at  the  free 
surface. 

vii.  Update  the  VOF  field: 

1.  For  each  fluid  element  calculate  the  element 
volume.  This  step  is  done  in  parallel  using 
the  list  of  fluid  elements  for  each  processor 
found  in  Step  l.g.ii. 

2.  For  each  fluid  element  find  the  apertures 
through  which  the  fluid  convects  to  each 
neighboring  element.  This  step  is  done  in 
parallel  using  the  list  of  fluid  elements  for 
each  processor  found  in  Step  l.g.ii. 

3.  For  each  fluid  element  use  the  apertures,  the 
element  volume,  the  element  current  VOF 
value,  and  the  element  nodal  velocities  to 
update  the  VOF  value  of  all  neighboring 
elements  by  finding  the  volume  of  fluid  that 
left  the  element  during  that  time  step  using 
Equations  24-28.  This  step  is  done  in  parallel 
using  the  list  of  fluid  elements  for  each 
processor  found  in  Step  l.g.ii.  Note  that  this 
step  depends  on  the  order  of  the  elements  in 
the  list  of  elements.  However,  since  the 
updates  of  the  VOF  field  between  solution 
time  steps  are  small,  therefore  this 
dependence  is  generally  very  small.  In  order 
to  assure  minimum  dependence  on  the 
elements’  order,  at  a  time  step  the  elements 
are  updated  from  first  to  last,  then  at  the  next 
time  step  they  are  updated  from  last  to  first. 


viii.Average  the  fluid  pressure  (This  step  eliminates 
the  pressure  checker-boarding  effect  and  allows 
use  of  equal  order  interpolation  for  both 
pressure  and  velocity). 

ix.Go  to  the  beginning  of  step  2. 

An  advantage  of  explicit  solution  procedures  is  that  they 
are  “embarrassingly”  parallel.  The  above  procedure 
achieves  near  linear  speed-up  with  the  number  of 
processors. 

7.  VALIDATION  STUDY 

A  test  fixture,  located  in  TARDEC’s  Simulation 
Laboratory  (TSL),  was  constructed  to  validate  the 
computational  code.  The  test  fixture  consists  of  an  oval 
tank  mounted  on  three  suspension  spring-dampers, 
which  are  in  turn  mounted  on  a  rigid  frame  (Figure  2). 
Connecting  rods  along  with  spherical  joints  are  used  to 
constrain  the  horizontal  (lateral  and  longitudinal)  motion 
of  the  springs,  such  that  they  move  nearly  vertically.  The 
frame  is  mounted  on  two  linear  hydraulic  actuators 
which  can  be  used  to  move  the  frame.  When  the 
actuators  move  the  same  way,  then  only  the  frame  pitch 
angle  is  changed.  When  the  actuators  move  in  an 
opposing  way,  then  the  roll  angle  of  the  frame  is 
changed.  The  test  fixture  is  designed  to  simulate  typical 
motions  that  a  ground  vehicle  is  subjected  to.  The  pitch 
motion  of  the  frame  simulates  both  sides  of  the  vehicle 
going  over  a  bump.  The  roll  motion  simulates  the  vehicle 
turning  or  going  over  a  bump  only  on  one  side  of  the 
vehicle. 


Figure  2  Test-fixture  model. 


The  tank  has  longitudinal  and  cross-section  baffles 
(Figure  3).  The  baffles  have  openings  near  the  bottom  of 
the  tank  to  equalize  the  liquid  level  when  the  tank  is  less 
than  half-full.  The  tank/baffles  geometry  and 
configuration  are  similar  to  typical  army  large-volume 
water/fuel  tanks. 


Figure  3  Tank  longitudinal  and  cross-sectional  baffles. 


The  linear  input  motion  of  the  two  actuators  is  measured 
using  LVDTs.  Also,  the  linear  motion  of  the  three 
suspension  springs  is  measured  using  LVDTs.  The 
signal  of  the  LVDTs  is  sampled  at  256  samples/sec.  Two 
cameras  are  mounted  on  the  ground  to  record  the 
motion  of  the  tank  and  sloshing  of  the  liquid  inside  the 
tank.  One  camera  shows  the  front  view  of  the  tank  and 
the  other  camera  shows  the  side  view  (Figure  4).  The 
cameras  are  set  to  capture  30  frames/sec. 


Figure  4  Front  (left)  and  side  (right)  views  of  the  tank  using  the  two 
ground  mounted  cameras. 

A  DIS  finite  element  model  of  the  test  fixture  is 

constructed.  The  model  consists  of  the  following 

components: 

•  An  oval  tank.  The  tank  dimensions  are:  0.324  m 
length,  width  0.428  m  and  height  0.276  m.  The  tank 
is  modeled  using  shell  elements. 

•  Three  linear  suspension  spring-dampers.  Each 
spring-damper  consist  of  a  linear  compression 
spring  in  parallel  with  a  linear  strut.  The  struts  were 
drained  of  fluid,  but  they  still  provide  a  small  amount 
of  viscous  damping  as  well  as  Coulomb  friction.  The 
physical  characteristics  of  the  springs-struts  are 
shown  in  Table  2.  The  struts  are  modeled  using 
cylindrical  joints. 

•  The  frame  is  modeled  as  a  rigid  body.  The  springs 
are  mounted  to  the  frame  using  spherical  joints. 

•  Two  linear  actuators  connected  to  the  frame. 

•  A  rigid  grounded  base. 

•  Connecting  rods.  Provide  horizontal  (lateral  and 
longitudinal)  stability  for  the  tank.  The  rods  are 
connected  to  the  tank  and  frame  using  spherical 
joints. 

•  17  spherical  joints. 


•  5  cylindrical  joints  located  at  the  3  suspension 

system  springs  and  the  two  actuators. 


Table  2  Physical  characteristics  of  the  test  fixture  spring-dampers. 


Left/Right  spring 

Middle  spring 

Stiffness 

3200  N/m 

8100  N/m 

Damping 

14  N  sec/m 

14  N  sec/m 

Friction  force 

6  N 

6  N 

The  test  fixture  is  simulated  with  an  empty  tank  and  with 
the  tank  half-full  with  water  (p=  1000  Kg/m3, 
p=  0.001  Kg/(m.  sec).  The  water  is  modeled  as 
incompressible  using  the  artificial  compressibility 
technique  with  an  artificial  sound  speed  factor  of  0.1  (i.e. 
the  artificial  sound  speed  in  the  water  is  taken  as 
1483m/sec  x  0.1  =  148.3  m/sec).  Due  to  the  use  of  large 
elements  near  the  solid  surface,  full  slip  boundary 
condition  at  the  wall  is  used  (s  =  1  in  Equation  23).  Thus, 
the  viscous  wall  friction  effects  are  assumed  to  be 
negligible.  Gravity  is  modeled  with  the  gravitational 
acceleration  taken  to  be  9.8  m/sec2  in  the  vertical 
direction. 


Table  3  Empty  tank  with  baffles  (Total  of  12  experiments). 


Pitch 

Roll 

Frequency 

(Hz) 

Amplitude 

(mm) 

Frequency 

(Hz) 

Amplitude 

(mm) 

Harmonic 

excitation 

1.5 

12 

2 

30 

2 

12 

3 

20 

3 

6 

4 

8 

Amplitude  (mm) 

Amplitude  (mm) 

Ramp 

excitation 

16 

10 

32 

20 

52 

52 

Table  4  Half-full  tank  with  baffles  and  half-full  tank  without  baffles 
(Total  of  24x2=48  experiments). 


Pitch 

Roll 

Stir 

Freq. 

(Hz) 

Amplitude 

(mm) 

Freq. 

(Hz) 

Amplitude 

(mm) 

Freq. 

(Hz) 

Amplitude 

(mm) 

Harmonic 

excitation 

1.5 

8 

1.5 

15 

1.5 

8 

1.5 

12 

1.5 

30 

1.5 

16 

2 

8 

2 

15 

2 

8 

2 

12 

2 

30 

2 

16 

3 

6 

3 

10 

3 

8 

3 

8 

3 

20 

4 

8 

Ramp 

excitation 

Amplitude  (mm) 

Amplitude  (mm) 

Amplitude  (mm) 

16 

10 

8/16 

32 

20 

16/32  mm 

A  test  matrix  consisting  of  60  experiments  was 
performed.  Tables  3  and  4  show  the  input  motion  of  the 
actuators  for  each  experiment.  The  input  motions  types 
include  pitch,  roll  and  stir  (combination  pitch  and  roll) 
excitations.  Recall  that  pitch  motion  means  that  the  two 
actuators  are  moving  in  the  same  way  (i.e.  in  phase)  and 
roll  excitation  means  that  the  two  actuators  are  moving 
in  an  opposing  way  (i.e.  180°  out  of  phase  with  each 
other  when  a  harmonic  excitation  is  used).  For  each 
motion  type  either  a  harmonic  excitation  or  a  ramp 
excitation  was  used.  For  the  harmonic  excitation  the 
frequency  and  amplitude  of  the  excitation  is  varied.  For 
the  ramp  excitation  just  the  amplitude  of  the  excitation  is 
varied. 


Figure  5  Comparison  between  experiment  and  DIS  simulation  for  an 
empty  tank  with  a  pitch  3  Hz,  6  mm  harmonic  excitation.  Top  2  graphs 
show  the  input  motion  of  the  actuator.  Bottom  three  graphs  show  the 
resulting  motion  of  the  3  suspension  springs. 


Figure  6  Comparison  between  experiment  and  DIS  simulation  for 
empty  tank  pitch  52  mm  ramp  excitation. 


Figure  7  Comparison  between  experiment  and  DIS  simulation  for 
empty  tank  roll  3  Hz,  20  mm  harmonic  roll  excitation. 


Figure  8  Comparison  between  experiment  and  DIS  simulation  for  half¬ 
full  tank  with  baffles  pitch  1.5  Hz,  12  mm  harmonic  excitation. 


Figure  9  Comparison  between  experiment  and  DIS  simulation  for  half¬ 
full  tank  with  baffles  pitch  2  Hz,  12  mm  harmonic  excitation. 


Figure  10  Comparison  between  experiment  and  DIS  simulation  for 
half-full  tank  with  baffles  pitch  3  Hz,  6  mm  harmonic  excitation. 


half-full  tank  with  baffles  pitch  3  Hz,  8  mm  harmonic  excitation.  Fi9ure  12  Comparison  between  experiment  and  DIS  simulation  for 

half-full  tank  with  baffles  roll  3  Hz,  20  mm  harmonic  excitation. 
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Figure  13  DIS  simulation  snapshots  of  the  front  view  of  the  half-filled  tank  with  baffles  at  3  Hz,  8  mm  pitch  excitation. 


Figure  14  Experiment  snapshots  of  the  front  view  of  the  half-filled  tank  with  baffles  at  3  Hz,  8  mm  pitch  excitation  taken  at  approximately  the  same  times 

as  the  snapshots  in  Figure  13. 


Figure  15  Simulation  and  experiment  snapshots  for  half-full  tank  with 
baffles  roll  2  Hz,  30  mm  harmonic  excitation. 


Table  3  shows  the  experiments  carried  out  with  an 
empty  tank.  Those  experiments  are  used  to  characterize 
the  suspension  spring  stiffness,  damping  and  friction 
properties,  as  well  as  to  validate  the  solid  multibody 
dynamics  solution.  Twelve  empty  tank  experiments  are 
carried  out.  Figure  5  shows  the  comparison  of  the 
experiment  and  simulation  results  for  a  3  Hz,  6  mm  pitch 
harmonic  excitation  with  an  empty  tank.  Since  both  the 
right  and  left  actuators  are  moving  the  same  amount, 
then  the  left  and  right  springs  should  also  move  the 
same  amount.  The  experiment  shows  that  the 
magnitude  of  motion  of  the  left  spring  is  about  10  mm 
while  the  magnitude  of  motion  of  the  right  spring  is  8 
mm.  The  DIS  simulation  predicts  that  the  magnitude  of 
motion  of  both  springs  is  about  7  mm.  The  experiment 
response  and  simulation  response  of  the  middle-spring 
are  practically  coincident. 

Figure  6  shows  a  comparison  of  the  experiment  and 
simulation  results  for  a  52  mm  pitch  ramp  excitation  of 
the  empty  tank.  Figure  7  shows  a  comparison  of  the 
experiment  and  simulation  results  for  a  3  Hz,  20  mm 
harmonic  roll  excitation  for  the  empty  tank.  Figures  5-7 
show  that  the  multibody  dynamics  model  of  the  test 
fixture  without  fluid  in  the  tank  can  predict  with  an 
average  difference  of  about  15%  the  response  of  the 
actual  test  fixture.  The  difference  between  the  results 
can  be  mostly  attributed  to: 

•  Non-linear  behavior  of  the  suspension  struts, 
including  non-linear  damping/friction  and 
clearances. 

•  Imprecision  of  the  spherical  joints  of  the  test  fixture. 
This  includes  joint  clearances  and  friction. 

This  is  evident  from  the  fact  that  in  Figure  5-7  the 
amplitude  of  motion  of  the  right  and  left  springs  are 
different  by  about  an  average  of  15%,  while  if  the  test 
fixture  joints  and  suspension  springs  were  ideal,  the 
difference  between  the  response  of  the  left  and  right 
springs  should  be  much  smaller  (~5%  similar  to  the 
difference  exhibited  in  the  simulation). 

Figures  8-1 1  show  a  comparison  of  the  experiment  and 
simulation  results  for  various  harmonic  pitch  excitations 


of  a  half-full  tank.  In  Figure  8  the  excitation  is  1 .5  Hz  and 
12  mm.  In  Figure  9  it  is  2  Hz,  12  mm.  In  Figure  10  it  is  3 
Hz,  6  mm  and  in  Figure  11  it  is  3  Hz,  8  mm.  Figure  12 
shows  a  comparison  of  the  experiment  and  simulation 
results  for  a  3  Hz,  20  mm  harmonic  roll  excitation  for  the 
half-full  tank.  Figures  8-12  show  that  the  DIS  simulation 
predicts  the  response  of  the  test  fixture  with  half-full  tank 
within  an  average  of  difference  of  15%  from  the 
experiments,  which  is  the  same  average  difference 
between  model  and  experiment  as  with  the  empty  tank 
runs.  This  suggests  that  the  difference  between  the 
experiment  and  simulation  is  mostly  due  to  the  non¬ 
linear  behavior  of  the  test  fixture  suspension  struts  and 
spherical  joints.  Figures  7  and  10  show  that  for  the  roll 
excitation  the  liquid  sloshing  in  the  tank  reduces  the 
motion  of  the  suspension  springs.  This  is  due  to  the  fact 
that  in  roll  motion,  the  fluid  acts  like  a  load  balancer,  thus 
limiting  the  transfer  of  forces  from  the  tank  to  the 
springs. 

Figure  13  shows  snapshots  of  the  front  view  of  the  half- 
filled  tank  with  baffles  at  3  Hz,  8  mm  pitch  excitation. 
Figure  14  shows  snapshots  taken  from  the  front  view 
experiment  camera  at  approximately  the  same  times  as 
the  snapshots  in  Figure  13.  Figures  13  and  14  show  that 
the  shape  and  mode  of  motion  of  the  liquid’s  free 
surface  predicted  by  the  model  is  approximately  the 
same  as  the  experiment.  Figures  11,  13  and  14  show 
that  at  3  Hz,  8  mm  harmonic  pitch  excitation,  the  liquid 
sloshing  undergoes  a  mode  change  from  a  straight  up- 
down  motion  along  with  the  tank  to  a  symmetric  sine 
sloshing  motion  along  the  tank  cross-section.  When  the 
liquid  undergoes  this  mode  change,  the  amplitude  of 
motion  of  the  suspension  springs  is  reduced  by  about 
25%.  Note  both  the  simulation  and  experiment  show  that 
this  mode  change  does  not  occur  at  3  Hz,  6mm 
harmonic  pitch  excitation  (Figure  10). 

Figure  15  shows  two  snapshots  of  the  front  view  of  the 
half-filled  tank  with  baffles  at  2  Hz,  30  mm  roll  excitation. 

8.  CONCLUDING  REMARKS 

A  finite  element  model  for  predicting  the  fully  coupled 
dynamic  response  of  flexible  multibody  systems  and 
liquid  sloshing  in  containers  was  presented.  The  model 
has  the  following  characteristics: 

•  Parallel  explicit  time-integration  solver. 

•  Library  of  accurate  large  rotation  finite  elements 
including:  truss,  beam,  shell  and  solid  elements.  The 
elements  only  use  Cartesian  coordinates  as  DOFs. 

•  The  fluid  mesh  is  modeled  using  a  very  light  and 
compliant  solid  mesh  which  allows  the  fluid  mesh  to 
move/deform  along  with  the  tank  using  the  Arbitrary 
Lagrangian-Eulerian  formulation. 

•  Acceptor-donor  VOF  algorithm  for  modeling  the 
fluid's  free-surface. 

•  The  motion  of  the  solid  and  fluid  is  referred  to  a 
global  inertial  Cartesian  reference  frame. 


•  A  total  Lagrangian  deformation  description  is  used 
for  the  solid  elements. 

•  The  penalty  technique  is  used  to  model  the  joints. 

A  validation  study  of  the  finite  element  model  was 
carried  out  using  a  specially  designed  test-fixture.  The 
study  shows  that  the  model  can  predict  reasonably  well 
(within  15%  on  average)  the  response  measured  on  the 
physical  test  fixture. 
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